/* ======Trend Inflation and the Nature of Structural Breaks in the New Keynesian Phillips Curve ======

This program estimates the unobserved components model for the NKPC with stochastic trend inflation
with two structural breaks. The transitory component z_t is constrained to be serially uncorrelated
and the variability of its shocks are constrained to be constant.

Written by: Pym Manopimoke

Model:

Pi_t = Pi_t^* + k_{S_t} Sum_{j=0}^{inf} beta^j E_{t-1}(x_{t+j}) + z_t
Pi_t^* = Pi_{t-1}^* + v_t, v_t|S_t ~iid.N(0,sig^2_{v,S_t))
x_t = phi_1 x_{t-1} + phi_2 x_{t-2} + u_t, u_t|D_t~iid.N(0,sig^2_{u,D_t})
z_t = epsilon_t, epsilon_t ~iid.N(0,sig^2_{epsilon})

u_t and epsilon_t are correlated with correlation coefficient rho.

S_t = {0,1,2}
Prob[ S(t)=0 | S(t-1)=0] = p00
Prob[ S(t)=1 | S(t-1)=1] = p11
Prob[ S(t)=2 | S(t-1)=2] = 1

sig^2_{u,D_t} = sig^2_{u,0}(1-D_t) + sig^2_{u,1}D_t
D_t = {0 if t<= 1984Q3, 1 otherwise}

========================================================================================================*/

new;
library optmum,PGRAPH;
pqgwin auto;
cls;
 
load data[250,11]=data.prn;  
data=data[2:250,.];
t=rows(data);
@Set the sample to 1952:1-2007:3@
startsmpl=21;
endsmpl =243; 

qtr= data[startsmpl:endsmpl,1];  
/*Column 2 for GDP chained-weighted inflation, 7 for CPI inflation, 8 for PCE inflation, 
9 for GDP Deflator*/  
inf=data[startsmpl:endsmpl,8];  @PCE inflation@
xx=data[startsmpl:endsmpl,5];  @CBO output gap@
t=rows(inf~xx);

@Define dummy variable series D_t to account for the Great Moderation@
GM = 130; @1984:2@
dmat = ones(GM,1)|zeros(t-GM,1); 

@================= Initialize Global Variables============@

tvp_no=4; 
start= 11;
t1=75;t2=125;

@Starting Values for PCE inflation@
prmtr_in=  4.3378770|        
           3.6575769|
           1.0308577|
           0.45933212|
           0.37570899| 
           1|
           0.29108135|
           0.71844084|
           -0.066901577|
           1.4701268|
           -3|
           0.022528669|
           0.0065056336|
           0.013872707;

 /*Robustness checks*/
 @Starting Values for CPI inflation@
 /*prmtr_in=   4.3178313|
              4.0997222|
              1.0309045| 
              0.45942266| 
              -0.47234109| 
              1.6868358 |
              0.26207587 |
              1.0131122 |
              -0.10269576| 
              1.4581428 |
              -4.3703318 |
              0.016971386 |
              0.026360364 |
              0.036512740; */
 
  @Starting values for GDP chain-weighted inflation@
  /*prmtr_in=  4.3634125  |
       4.0198254  |
       1.0328900  |
      0.45902465  |
     -0.38108505  |
       1.3153017  |
      0.26964884  |
      0.61142235  |
     -0.33543943  |
       1.4196333  |
      -4.6105111  |
     0.015869413  |
    -0.013394352  |
     0.012067873  ;  */
  
{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);
PRM_FNL=TRANS(xout);           
xout;
hout0=hessp(&lik_fcn,xout);
hout0;
hout=inv(hout0);
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl));        @Standard errors of the estimated coefficients@
 
{B_TT,PR_TT,PR_TL,stdprederr,xfor,confint} =FILTER(XOUT);
smooth3=SMOOTH(pr_tt,pr_tl,prm_fnl[1:2,1]); @Pr[S_t|Y_T], Smoothed probabilities@

"== FINAL OUTPUT ===============================";
"Likelihood value is ";; -fout;
"Estimated parameters (s.e.) tstat are:"; 
prm_fnl~sd_fnl~prm_fnl./sd_fnl;
"===============================================";
 
_pdate=0;
_pltype = {4 6 3};
_plegctl={4 4 6 5};
_plegstr="Regime1\000Regime2\000Regime3";
title("Smoothed probabilities of regimes");
XY(qtr[start:t],SMOOTH3); 
 
_pltype = {6 3 4 4};
_plegctl={2 4 5.5 4.8};
_plegstr="Trend Inflation\000Actual Inflation\00090\% Confidence Bands";
title("Actual inflation and trend inflation");
xy(qtr[start:t],b_tt[.,1]~inf[start:t]~confint);
 
_pltype = {6 3};
_plegctl={2 4 6 5};
_plegstr="Inflation Cycle\000Inflation Gap";
title("Inflation cycle and inflation gap");
XY(qtr[start:t],xfor~inf[start:t]-b_tt[.,1]);
    
output file=UC1.out reset;
output off;
end;

@=================================================================================@

PROC LIK_FCN(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2, 
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22;  
 
EXTERNAL PROC TRANS, V_PROB;
 
    PRMTR=TRANS(PRMTR1);
        
    P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
    P11=PRMTR[2,1];    @Pr[St=1/St-1=1]@
    P22=1;             @Pr[St=2/St-1=2]@
                
    siguGM0 = prmtr[3,1];
    siguGM1 = prmtr[4,1];
    sigv= prmtr[5:7,1]; 
    sigep = prmtr[8,1]; 
    rho = prmtr[9,1];
                    
    phi1 = prmtr[10,1];
    phi2 = prmtr[11,1];
    psi=0;    
    k=prmtr[12:14,1];
    
    A=1~0;      @selection matrix@
    H = 1~1~0~0|0~0~1~0;
               
    B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
    P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
    B_LL1=inf[t1]|0|xx[1]|0; 
    P_LL1 = eye(4)*10000;       
    B_LL2=inf[t2]|0|xx[1]|0; 
    P_LL2 = eye(4)*10000;      
                
    subF = phi1~phi2| 1~0; 
    subQ= siguGM0^2~0|0~0;
    vecp = inv(eye(4)-subF.*.subf)*vec(subq);
    vv = reshape(vecp',2,2);
    p_ll0[3:4,3:4] = vv;
    p_ll1[3:4,3:4] = vv;
    p_ll2[3:4,3:4] = vv;

    p_ll0[2,2] = sigep^2/(1-psi^2);
    p_ll1[2,2] = p_ll0[2,2];
    p_ll2[2,2] = p_ll0[2,2];

    PROB__0 = 1; @Unconditional probability for St=0@
    PROB__1 = 0; @Unconditional probability for St=1@
    PROB__2 = 0; @Unconditional probability for St=2@
 
/*  START ITERATION */
/*===================================================================== */

    LIKV=0.0;
    J_ITER = 1;
    do until J_ITER>T;

        yt=inf[j_iter]|xx[j_iter];
        dmt = dmat[j_iter];
        sigu=siguGM0*dmt + siguGM1*(1-dmt);
 
        Q0=    sigv[1]^2~0~0~0|
               0~sigep^2~rho*sigu*sigep~0|
               0~rho*sigu*sigep~sigu^2~0|
               0~0~0~0;
        Q1=    sigv[2]^2~0~0~0|
               0~sigep^2~rho*sigu*sigep~0|
               0~rho*sigu*sigep~sigu^2~0|
               0~0~0~0;
        Q2=    sigv[3]^2~0~0~0|
               0~sigep^2~rho*sigu*sigep~0|
               0~rho*sigu*sigep~sigu^2~0|
               0~0~0~0;
        F0 =   1~0~0~0 |
               0~psi~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F1 =   1~0~0~0 |
               0~psi~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F2 =   1~0~0~0 |
               0~psi~0~0|
               0~0~phi1~phi2|
               0~0~1~0;

        FF= phi1~phi2|1~0;
        xfor = k[1]*A*(inv(eye(2)-FF));  @infinite sum term in cyclical component@        
        ome0 =0~0~xfor| 0~0~0~0;
        xfor = k[2] *A*(inv(eye(2)-FF)); 
        ome1 =0~0~xfor| 0~0~0~0;
        xfor = k[3] *A*(inv(eye(2)-FF)); 
        ome2 =0~0~xfor| 0~0~0~0;

        /* KALMAN FILTER */
        @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
        B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
        B_TL01 = F1 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
        B_TL02 = F2 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=2@
        B_TL10 = F0 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
        B_TL11 = F1 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@
        B_TL12 = F2 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=2@
        B_TL20 = F0 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=0@
        B_TL21 = F1 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=1@
        B_TL22 = F2 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=2@
        
        @ P_tl(i,j),  S(t-1)=i, S(t)=j@
        P_TL00 = F0 * P_LL0 * F0' + Q0;
        P_TL01 = F1 * P_LL0 * F1' + Q1;
        P_TL02 = F2 * P_LL0 * F2' + Q2;
        P_TL10 = F0 * P_LL1 * F0' + Q0;
        P_TL11 = F1 * P_LL1 * F1' + Q1;
        P_TL12 = F2 * P_LL1 * F2' + Q2;
        P_TL20 = F0 * P_LL2 * F0' + Q0;
        P_TL21 = F1 * P_LL2 * F1' + Q1;
        P_TL22 = F2 * P_LL2 * F2' + Q2;
 
        @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
        F_CAST00= yt- H * B_TL00 - ome0* B_TL00;
        F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
        F_CAST02= yt- H * B_TL02 - ome2* B_TL02;
        F_CAST10= yt- H * B_TL10 - ome0* B_TL10;
        F_CAST11= yt- H * B_TL11 - ome1* B_TL11;
        F_CAST12= yt- H * B_TL12 - ome2* B_TL12;
        F_CAST20= yt- H * B_TL20 - ome0* B_TL20;
        F_CAST21= yt- H * B_TL21 - ome1* B_TL21;
        F_CAST22= yt- H * B_TL22 - ome2* B_TL22;
 
        @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
        SS00= H * P_TL00 * H' ; 
        SS01= H*  P_TL01 * H' ;
        SS02= H*  P_TL02 * H' ;
        SS10= H * P_TL10 * H' ;
        SS11= H * P_TL11 * H' ;
        SS12= H * P_TL12 * H' ;
        SS20= H * P_TL20 * H' ;
        SS21= H * P_TL21 * H' ;
        SS22= H * P_TL22 * H' ;
 
        /*Updating*/
        B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
        B_TT01 = B_TL01 + (P_TL01 * H' *inv(SS01)) * F_CAST01;
        B_TT02 = B_TL02 + (P_TL02 * H' *inv(SS02)) * F_CAST02;
        B_TT10 = B_TL10 + (P_TL10 * H' *inv(SS10)) * F_CAST10;
        B_TT11 = B_TL11 + (P_TL11 * H' *inv(SS11)) * F_CAST11;
        B_TT12 = B_TL12 + (P_TL12 * H' *inv(SS12)) * F_CAST12;
        B_TT20 = B_TL20 + (P_TL20 * H' *inv(SS20)) * F_CAST20;
        B_TT21 = B_TL21 + (P_TL21 * H' *inv(SS21)) * F_CAST21;
        B_TT22 = B_TL22 + (P_TL22 * H' *inv(SS22)) * F_CAST22;

        P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
        P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
        P_TT02 = (EYE(tvp_no) - (P_TL02 * H' *inv(SS02) * H )) * P_TL02;
        P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
        P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
        P_TT12 = (EYE(tvp_no) - (P_TL12 * H' *inv(SS12) * H )) * P_TL12;
        P_TT20 = (EYE(tvp_no) - (P_TL20 * H' *inv(SS20) * H )) * P_TL20;
        P_TT21 = (EYE(tvp_no) - (P_TL21 * H' *inv(SS21) * H )) * P_TL21;
        P_TT22 = (EYE(tvp_no) - (P_TL22 * H' *inv(SS22) * H )) * P_TL22;
       
        /* HAMILTON FILTER */   

        PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
        PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
        PR_VL02=0;
        PR_VL10=0;
        PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
        PR_VL12=V_PROB(F_CAST12,SS12)*(1-p11)*PROB__1;
        PR_VL20=0;
        PR_VL21=0;
        PR_VL22=V_PROB(F_CAST22,SS22)* p22 *PROB__2;
        
        @CONDITIONAL DENSITY TIMES WEIGHT@

        PR_VAL=PR_VL00+PR_VL01+PR_VL02+PR_VL10+PR_VL11+PR_VL12+PR_VL20+PR_VL21+PR_VL22;
                                  
        @Pr[St,St-1/Yt]@
        PRO_00=PR_VL00/PR_VAL;       
        PRO_01=PR_VL01/PR_VAL;
        PRO_02=PR_VL02/PR_VAL;
        PRO_10=PR_VL10/PR_VAL;
        PRO_11=PR_VL11/PR_VAL;
        PRO_12=PR_VL12/PR_VAL;
        PRO_20=PR_VL20/PR_VAL;
        PRO_21=PR_VL21/PR_VAL;
        PRO_22=PR_VL22/PR_VAL;
        PROB__0=PRO_00+PRO_10+PRO_20;       @Pr[St=0/Yt]@
        PROB__1=PRO_01+PRO_11+PRO_21;       @Pr[St=1/Yt]@
        PROB__2=PRO_02+PRO_12+PRO_22;       @Pr[St=2/Yt]@

       /* COLLAPSING TERMS */
      
       if(prob__0!=0);
       B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 + PRO_20*B_TT20)/PROB__0;
       endif;
       if(prob__1!=0);
       B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 + PRO_21*B_TT21)/PROB__1;
       endif;
       if(prob__2!=0);
       B_LL2=(PRO_02*B_TT02 + PRO_12*B_TT12 + PRO_22*B_TT22)/PROB__2;
       endif;
       if(prob__0!=0);
       P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
       PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)')+
       PRO_20*(P_TT20+(B_LL0-B_TT20)*(B_LL0-B_TT20)'))/PROB__0;
       endif;
       if(prob__1!=0);
       P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
       PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)')+
       PRO_21*(P_TT21+(B_LL1-B_TT21)*(B_LL1-B_TT21)'))/PROB__1;
       endif;
       if(prob__2!=0);
       P_LL2=(PRO_02*(P_TT02+(B_LL2-B_TT02)*(B_LL2-B_TT02)')+
       PRO_12*(P_TT12+(B_LL2-B_TT12)*(B_LL2-B_TT12)')+
       PRO_22*(P_TT22+(B_LL2-B_TT22)*(B_LL2-B_TT22)'))/PROB__2;
       endif;
       if J_ITER < START; goto skip; endif;
       LIKV = LIKV -LN(PR_VAL);
       skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
  
RETP(LIKV);
ENDP;
 
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc TRANS(prmtr1); @constraining values of reg. coeff.@
local prm_new,d0,d1;

prm_new=prmtr1;  
 
prm_new[1:2,1]=EXP(prmtr1[1:2,1]) ./ (1 + EXP(prmtr1[1:2,1])); 
prm_new[3:8,1] = sqrt(prmtr1[3:8,1]^2);
prm_new[9] =prmtr1[9,1]./(1+abs(prmtr1[9,1]));
d0 =prmtr1[10,1]./(1+abs(prmtr1[10,1]));
d1 =(1-abs(d0))*prmtr1[11,1]./(1+abs(prmtr1[11,1]))+ abs(d0)-d0^2;
prm_new[10,1]=2*d0;
prm_new[11,1]=-1*(d0^2 + d1) ;

retp(prm_new);
endp;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc V_PROB(EV,HE);  @ CALCULATES Pr[Yt/St,Yt-1]@
local VAL;

VAL=(1/SQRT(2*PI*det(HE)))*EXP(-0.5*EV'inv(HE)*EV);

RETP(VAL);
ENDP;
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@

proc (6) = FILTER(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor0,xfor1,xfor2,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,confint,lowerCI,upperCI,stderr,stdprederr,xforecast,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2,B_M,PR_TTM,PR_TLM,
pr_tl00,pr_tl01, pr_tl02, pr_tl10,pr_tl11,pr_tl12,pr_tl20,pr_tl21,pr_tl22,
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22;  

EXTERNAL PROC TRANS, V_PROB;
    
    PRMTR=TRANS(PRMTR1);
 
    P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
    P11=PRMTR[2,1];    @Pr[St=1/St-1=1]@
    P22=1;             @Pr[St=2/St-1=2]@
                
    siguGM0 = prmtr[3,1];
    siguGM1 = prmtr[4,1];
    sigv= prmtr[5:7,1]; 
    sigep = prmtr[8,1]; 
    rho = prmtr[9,1];
                    
    phi1 = prmtr[10,1];
    phi2 = prmtr[11,1];
    psi=0;    
    k=prmtr[12:14,1];
    
    A=1~0;      @selection matrix@
    H = 1~1~0~0|0~0~1~0;
               
    B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
    P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
    B_LL1=inf[t1]|0|xx[1]|0; 
    P_LL1 = eye(4)*10000;       
    B_LL2=inf[t2]|0|xx[1]|0; 
    P_LL2 = eye(4)*10000;      
                
    subF = phi1~phi2| 1~0; 
    subQ= siguGM0^2~0|0~0;
    vecp = inv(eye(4)-subF.*.subf)*vec(subq);
    vv = reshape(vecp',2,2);
    p_ll0[3:4,3:4] = vv;
    p_ll1[3:4,3:4] = vv;
    p_ll2[3:4,3:4] = vv;

    p_ll0[2,2] = sigep^2/(1-psi^2);
    p_ll1[2,2] = p_ll0[2,2];
    p_ll2[2,2] = p_ll0[2,2];

    PROB__0 = 1; @Unconditional probability for St=0@
    PROB__1 = 0; @Unconditional probability for St=1@
    PROB__2 = 0; @Unconditional probability for St=2@
 
    B_M=ZEROS(T,tvp_no);  @will store (1,1)element of B_TT@
    PR_TTM=ZEROS(T,3);  @will store Pr[S_t=0|Y_t]@
    PR_TLM=ZEROS(T,3);  @will store Pr[S_t=0|Y_{t-1}]@
    stdprederr=zeros(t,1); @standardized prediction error for diagnostic tests@
    confint=zeros(t,2);  @90% confidence interval bands for trend inflation estimates@ 
    xforecast=zeros(t,1); @stores the infinite sum term in the cyclical component@

/*  START ITERATION */
/*===================================================================== */

   LIKV=0.0;
   J_ITER = 1;
   do until J_ITER>T;
        
        yt=inf[j_iter]|xx[j_iter];
        dmt = dmat[j_iter];
        sigu=siguGM0*dmt + siguGM1*(1-dmt);
           
        pr_tl00 = PROB__0*p00;         
        pr_tl01 = PROB__0*(1-p00);
        pr_tl02 = 0;
        pr_tl10 = 0;
        pr_tl11 = PROB__1*p11;
        pr_tl12 = PROB__1*(1-p11);
        pr_tl20 = 0;
        pr_tl21 = 0;
        pr_tl22 = p22 *PROB__2;
              
        PR_TLM[J_ITER,.]= (pr_tl00 +pr_tl10+pr_tl20)
                         ~(pr_tl01 +pr_tl11+pr_tl21)
                         ~(pr_tl02 +pr_tl12+pr_tl22); 
 
        Q0=    sigv[1]^2~0~0~0|
               0~sigep^2~rho*sigu*sigep~0|
               0~rho*sigu*sigep~sigu^2~0|
               0~0~0~0;
        Q1=    sigv[2]^2~0~0~0|
               0~sigep^2~rho*sigu*sigep~0|
               0~rho*sigu*sigep~sigu^2~0|
               0~0~0~0;
        Q2=    sigv[3]^2~0~0~0|
               0~sigep^2~rho*sigu*sigep~0|
               0~rho*sigu*sigep~sigu^2~0|
               0~0~0~0;
        F0 =   1~0~0~0 |
               0~psi~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F1 =   1~0~0~0 |
               0~psi~0~0|
               0~0~phi1~phi2|
               0~0~1~0;
        F2 =   1~0~0~0 |
               0~psi~0~0|
               0~0~phi1~phi2|
               0~0~1~0;

        FF= phi1~phi2|1~0;
        xfor0 = k[1]*A*(inv(eye(2)-FF));  @infinite sum term in cyclical component@        
        ome0 =0~0~xfor0| 0~0~0~0;
        xfor1 = k[2] *A*(inv(eye(2)-FF)); 
        ome1 =0~0~xfor1| 0~0~0~0;
        xfor2 = k[3] *A*(inv(eye(2)-FF)); 
        ome2 =0~0~xfor2| 0~0~0~0;

        /* KALMAN FILTER */
        @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
        B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
        B_TL01 = F1 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
        B_TL02 = F2 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=2@
        B_TL10 = F0 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
        B_TL11 = F1 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@
        B_TL12 = F2 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=2@
        B_TL20 = F0 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=0@
        B_TL21 = F1 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=1@
        B_TL22 = F2 * B_LL2;  @WHEN S_{t-1}=2, S_{t}=2@
        
        @ P_tl(i,j),  S(t-1)=i, S(t)=j@
        P_TL00 = F0 * P_LL0 * F0' + Q0;
        P_TL01 = F1 * P_LL0 * F1' + Q1;
        P_TL02 = F2 * P_LL0 * F2' + Q2;
        P_TL10 = F0 * P_LL1 * F0' + Q0;
        P_TL11 = F1 * P_LL1 * F1' + Q1;
        P_TL12 = F2 * P_LL1 * F2' + Q2;
        P_TL20 = F0 * P_LL2 * F0' + Q0;
        P_TL21 = F1 * P_LL2 * F1' + Q1;
        P_TL22 = F2 * P_LL2 * F2' + Q2;
 
        @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
        F_CAST00= yt- H * B_TL00 - ome0* B_TL00;
        F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
        F_CAST02= yt- H * B_TL02 - ome2* B_TL02;
        F_CAST10= yt- H * B_TL10 - ome0* B_TL10;
        F_CAST11= yt- H * B_TL11 - ome1* B_TL11;
        F_CAST12= yt- H * B_TL12 - ome2* B_TL12;
        F_CAST20= yt- H * B_TL20 - ome0* B_TL20;
        F_CAST21= yt- H * B_TL21 - ome1* B_TL21;
        F_CAST22= yt- H * B_TL22 - ome2* B_TL22;
 
        @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
        SS00= H * P_TL00 * H' ; 
        SS01= H*  P_TL01 * H' ;
        SS02= H*  P_TL02 * H' ;
        SS10= H * P_TL10 * H' ;
        SS11= H * P_TL11 * H' ;
        SS12= H * P_TL12 * H' ;
        SS20= H * P_TL20 * H' ;
        SS21= H * P_TL21 * H' ;
        SS22= H * P_TL22 * H' ;
 
        /*Updating*/
        B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
        B_TT01 = B_TL01 + (P_TL01 * H' *inv(SS01)) * F_CAST01;
        B_TT02 = B_TL02 + (P_TL02 * H' *inv(SS02)) * F_CAST02;
        B_TT10 = B_TL10 + (P_TL10 * H' *inv(SS10)) * F_CAST10;
        B_TT11 = B_TL11 + (P_TL11 * H' *inv(SS11)) * F_CAST11;
        B_TT12 = B_TL12 + (P_TL12 * H' *inv(SS12)) * F_CAST12;
        B_TT20 = B_TL20 + (P_TL20 * H' *inv(SS20)) * F_CAST20;
        B_TT21 = B_TL21 + (P_TL21 * H' *inv(SS21)) * F_CAST21;
        B_TT22 = B_TL22 + (P_TL22 * H' *inv(SS22)) * F_CAST22;

        P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
        P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
        P_TT02 = (EYE(tvp_no) - (P_TL02 * H' *inv(SS02) * H )) * P_TL02;
        P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
        P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
        P_TT12 = (EYE(tvp_no) - (P_TL12 * H' *inv(SS12) * H )) * P_TL12;
        P_TT20 = (EYE(tvp_no) - (P_TL20 * H' *inv(SS20) * H )) * P_TL20;
        P_TT21 = (EYE(tvp_no) - (P_TL21 * H' *inv(SS21) * H )) * P_TL21;
        P_TT22 = (EYE(tvp_no) - (P_TL22 * H' *inv(SS22) * H )) * P_TL22;
           
       /* HAMILTON FILTER */   

        PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
        PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
        PR_VL02=0;
        PR_VL10=0;
        PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
        PR_VL12=V_PROB(F_CAST12,SS12)*(1-p11)*PROB__1;
        PR_VL20=0;
        PR_VL21=0;
        PR_VL22=V_PROB(F_CAST22,SS22)* p22 *PROB__2;
        
        @CONDITIONAL DENSITY TIMES WEIGHT@

        PR_VAL=PR_VL00+PR_VL01+PR_VL02+PR_VL10+PR_VL11+PR_VL12+PR_VL20+PR_VL21+PR_VL22;
                                  
        @Pr[St,St-1/Yt]@
        PRO_00=PR_VL00/PR_VAL;       
        PRO_01=PR_VL01/PR_VAL;
        PRO_02=PR_VL02/PR_VAL;
        PRO_10=PR_VL10/PR_VAL;
        PRO_11=PR_VL11/PR_VAL;
        PRO_12=PR_VL12/PR_VAL;
        PRO_20=PR_VL20/PR_VAL;
        PRO_21=PR_VL21/PR_VAL;
        PRO_22=PR_VL22/PR_VAL;
        
        PROB__0=PRO_00+PRO_10+PRO_20;       @Pr[St=0/Yt]@
        PROB__1=PRO_01+PRO_11+PRO_21;       @Pr[St=1/Yt]@
        PROB__2=PRO_02+PRO_12+PRO_22;       @Pr[St=2/Yt]@
        PR_TTM[J_ITER,.]=PROB__0~prob__1~prob__2;
       /* COLLAPSING TERMS */
      
       if(prob__0!=0);
       B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 + PRO_20*B_TT20)/PROB__0;
       endif;
       if(prob__1!=0);
       B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 + PRO_21*B_TT21)/PROB__1;
       endif;
       if(prob__2!=0);
       B_LL2=(PRO_02*B_TT02 + PRO_12*B_TT12 + PRO_22*B_TT22)/PROB__2;
       endif;
       if(prob__0!=0);
       P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
       PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)')+
       PRO_20*(P_TT20+(B_LL0-B_TT20)*(B_LL0-B_TT20)'))/PROB__0;
       endif;
       if(prob__1!=0);
       P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
       PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)')+
       PRO_21*(P_TT21+(B_LL1-B_TT21)*(B_LL1-B_TT21)'))/PROB__1;
       endif;
       if(prob__2!=0);
       P_LL2=(PRO_02*(P_TT02+(B_LL2-B_TT02)*(B_LL2-B_TT02)')+
       PRO_12*(P_TT12+(B_LL2-B_TT12)*(B_LL2-B_TT12)')+
       PRO_22*(P_TT22+(B_LL2-B_TT22)*(B_LL2-B_TT22)'))/PROB__2;
       endif;
 
       @Estimates of the unobserved components@
       B_TT=   B_TT00*PRO_00+B_TT01*PRO_01+B_TT02*PRO_02+
                B_TT10*PRO_10+B_TT11*PRO_11+B_TT12*PRO_12+
                B_TT20*PRO_20+B_TT21*PRO_21+B_TT22*PRO_22;
       B_M[J_ITER,.]=B_TT';

       @90% Confidence bands for trend inflation estimates@
       lowerCI=(B_TT00[1,1]-1.645*sqrt(P_TT00[1,1]))*PR_tl00+
               (B_TT01[1,1]-1.645*sqrt(P_TT01[1,1]))*PR_tl01+
               (B_TT02[1,1]-1.645*sqrt(P_TT02[1,1]))*PR_tl02+
               (B_TT10[1,1]-1.645*sqrt(P_TT10[1,1]))*PR_tl10+
               (B_TT11[1,1]-1.645*sqrt(P_TT11[1,1]))*PR_tl11+
               (B_TT12[1,1]-1.645*sqrt(P_TT12[1,1]))*PR_tl12+
               (B_TT20[1,1]-1.645*sqrt(P_TT20[1,1]))*PR_tl20+
               (B_TT21[1,1]-1.645*sqrt(P_TT21[1,1]))*PR_tl21+
               (B_TT22[1,1]-1.645*sqrt(P_TT22[1,1]))*PR_tl22;
       upperCI=(B_TT00[1,1]+1.645*sqrt(P_TT00[1,1]))*PR_tl00+
               (B_TT01[1,1]+1.645*sqrt(P_TT01[1,1]))*PR_tl01+
               (B_TT02[1,1]+1.645*sqrt(P_TT02[1,1]))*PR_tl02+
               (B_TT10[1,1]+1.645*sqrt(P_TT10[1,1]))*PR_tl10+
               (B_TT11[1,1]+1.645*sqrt(P_TT11[1,1]))*PR_tl11+
               (B_TT12[1,1]+1.645*sqrt(P_TT12[1,1]))*PR_tl12+
               (B_TT20[1,1]+1.645*sqrt(P_TT20[1,1]))*PR_tl20+
               (B_TT21[1,1]+1.645*sqrt(P_TT21[1,1]))*PR_tl21+
               (B_TT22[1,1]+1.645*sqrt(P_TT22[1,1]))*PR_tl22;
       confint[j_iter,.] = lowerCI~upperCI;

       @Standardized prediction error@
       stderr= (f_cast00[1,1]/sqrt(ss00[1,1]))*pr_tl00+ (f_cast01[1,1]/sqrt(ss01[1,1]))*pr_tl01+
               (f_cast02[1,1]/sqrt(ss02[1,1]))*pr_tl02+ (f_cast10[1,1]/sqrt(ss10[1,1]))*pr_tl10+
               (f_cast11[1,1]/sqrt(ss11[1,1]))*pr_tl11+ (f_cast12[1,1]/sqrt(ss12[1,1]))*pr_tl12+ 
               (f_cast20[1,1]/sqrt(ss20[1,1]))*pr_tl20+ (f_cast21[1,1]/sqrt(ss21[1,1]))*pr_tl21+
               (f_cast22[1,1]/sqrt(ss22[1,1]))*pr_tl22;
       STDPREDERR[j_iter,.] = stderr;
       
       @infinite sum term in the cyclical component for inflation@
       xforecast[j_iter,.]= (xfor0 *B_TL00[3:4,1])*pr_tl00 +
                            (xfor1 *B_TL01[3:4,1])*pr_tl01 +
                            (xfor2 *B_TL02[3:4,1])*pr_tl02 +
                            (xfor0 *B_TL10[3:4,1])*pr_tl10 +
                            (xfor1 *B_TL11[3:4,1])*pr_tl11 +
                            (xfor2 *B_TL12[3:4,1])*pr_tl12 + 
                            (xfor0 *B_TL20[3:4,1])*pr_tl20 +
                            (xfor1 *B_TL21[3:4,1])*pr_tl21 +
                            (xfor2 *B_TL22[3:4,1])*pr_tl22;

       if J_ITER < START; goto skip; endif;
       LIKV = LIKV -LN(PR_VAL);
       skip:
       J_ITER = J_ITER+1;
       endo;
  
/*  END ITERATION    */
/*===================*/

retp(B_M[START:T,.],PR_TTM[START:T,.],PR_TLM[START:T,.],stdprederr[start:t,.],xforecast[start:t,.],confint[start:t,.]);
ENDP;


@==========================================================@

@Kim's smoothing@
proc SMOOTH(prtt,prtl,prmtr0);

local   pr_sm1,pr_sm2,pr_sm3, j_iter, TT,prmtr, pr_tr, p11,p12,p13,p21,p22,p23,p31,p32,p33,
        pr_tt1,pr_tt2,pr_tt3,pr_tl1,pr_tl2,pr_tl3,pr_sm11,pr_sm12,pr_sm13,pr_sm21,pr_sm22,pr_sm23,pr_sm31,pr_sm32,pr_sm33;

      TT=ROWS(PRTT);
      prmtr=prmtr0;
      p11=prmtr[1,1]; 
      p12=1-p11; 
      p13=0; 
      p21=0; 
      p22=prmtr[2,1]; 
      p23=1-p22;
      p31=0; 
      p32=0; 
      p33=1;
 
      pr_tt1=prtt[.,1];
      pr_tt2=prtt[.,2];
      pr_tt3=prtt[.,3];
      pr_tl1=prtl[.,1];
      pr_tl2=prtl[.,2];
      pr_tl3=prtl[.,3];
           
      pr_sm1=pr_tt1;      @ pr_sm1 will contain Pr[S_t=1|Y_T]@
      pr_sm2=pr_tt2;      @ pr_sm2 will contain Pr[S_t=2|Y_T]@
      pr_sm3=pr_tt3;      @ pr_sm3 will contain Pr[S_t=3|Y_T]@
          

     j_iter=TT-1;
     do until j_iter < 1;
  
     if pr_tl1[j_iter+1,1]<0.000001; pr_tl1[j_iter+1,1]=0.0000001; endif;
     if pr_tl2[j_iter+1,1]<0.000001; pr_tl2[j_iter+1,1]=0.0000001; endif;
     if pr_tl3[j_iter+1,1]<0.000001; pr_tl3[j_iter+1,1]=0.0000001; endif;
 
        pr_sm11=pr_sm1[j_iter+1,1]*p11*  pr_tt1[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm12=pr_sm2[j_iter+1,1]*p12*  pr_tt1[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm13=pr_sm3[j_iter+1,1]*p13*  pr_tt1[j_iter,1]/ pr_tl3[j_iter+1,1];
 
        pr_sm21=pr_sm1[j_iter+1,1]*p21*  pr_tt2[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm22=pr_sm2[j_iter+1,1]*p22*  pr_tt2[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm23=pr_sm3[j_iter+1,1]*p23*  pr_tt2[j_iter,1]/ pr_tl3[j_iter+1,1];
      
        pr_sm31=pr_sm1[j_iter+1,1]*p31*  pr_tt3[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm32=pr_sm2[j_iter+1,1]*p32*  pr_tt3[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm33=pr_sm3[j_iter+1,1]*p33*  pr_tt3[j_iter,1]/ pr_tl3[j_iter+1,1];
       
        pr_sm1[j_iter,1]=pr_sm11+pr_sm12+pr_sm13;
        pr_sm2[j_iter,1]=pr_sm21+pr_sm22+pr_sm23;
        pr_sm3[j_iter,1]=pr_sm31+pr_sm32+pr_sm33;
    
     j_iter=j_iter -1;
     endo;


retp(pr_sm1~pr_sm2~pr_sm3);  
endp;

/*=================================================================================================================*/
 